#ifndef redlesiondetectionH
#define redlesiondetectionH

#include "Image.H"
#include "PixelIterator.H"
#include "Util.H"
#include "IO.H"
#include "graymorphology.h"
#include "funduspreprocessing.h"

#define _USE_MATH_DEFINES 
#include <math.h>


template<class T>
inline int round(const T& i) {
  return (int)(i + (T)0.5);
}

template <class Tin>
class RedLesionDetector{
        private:
        //builds a linear structuring element rotated to a certain
        //orientation (in degrees)
        void buildSE(pip::Image<unsigned char>& imSE, int orientation, int length);

        int nrOfOrientations;
        int elementLength;

        public:
        RedLesionDetector() : nrOfOrientations(6), elementLength(5){};
        //Black top-hat transform:
        // Do a closing with a linear structuring element at different
        //  rotational orientations take the minimum of the resulting images and
        //  subtract it from the input image
        int run(const pip::Image<Tin>& imIn, pip::Image<Tin>& imOut);
        int run(const pip::Image<pip::Color<unsigned char> >& imIn, pip::Image<Tin>& imOut);

        int setNrOfOrientations(int _nrOfOrientations){
                nrOfOrientations = _nrOfOrientations;
                return nrOfOrientations;
        }
        int setElementLength(int _elementLength){
                elementLength= _elementLength;
                return elementLength;
        }
};

//---------------------------------------------------------------------------
//////////
// Rotation around z-axis using nearest neighbour interpolation. only 2D
template <class Type>
int rotateNearestZ(pip::Image<Type>& imIn,
const float fRotAngleZ,
const Type theBackground) {

        // Initialization:
        int iWidth    = imIn.dim(0);
	    int iHeight   = imIn.dim(1);
		pip::Image<Type> imOut(imIn);
        int iMaxX = iWidth - 1;
        int iMaxY = iHeight - 1;
        const float fMaxX = (float)iMaxX;
        const float fMaxY = (float)iMaxY;

        // center of rotation
        const float fXc = fMaxX/2.f;
        const float fYc = fMaxY/2.f;

        // Precompute the sinus and cosines of the rotation angle. Notice
        // that for the eventual rotation loop, the inverse of the rotation
        // matrix is needed. For a rotation around the z-axis only, this is
        // the same as taking the negative angle:
        const float fAngleZ = fRotAngleZ*M_PI/180.0;
        const float fCosAZ = cos(-fAngleZ);
        const float fSinAZ = sin(-fAngleZ);

        if (imOut.dimension()!= imIn.dimension()) {
                imOut.resize(imIn.dimension());
        }

        // Rotation:
        for (int y=0; y<iHeight; ++y)
        for (int x=0; x<iWidth; ++x) {
                const float fdx = (float)x - fXc;
                const float fdy = (float)y - fYc;
                const int iX = (int)round(fXc + fdx*fCosAZ - fdy*fSinAZ);
                const int iY = (int)round(fYc + fdx*fSinAZ + fdy*fCosAZ);
                if ((iX < 0) || (iX > iMaxX) ||
                (iY < 0) || (iY > iMaxY))
                imOut(x,y) = theBackground;
                else
                imOut(x,y) = imIn(iX,iY);
        }
		imIn=imOut;

        return 1;
}

//---------------------------------------------------------------------------
//builds a linear structuring element rotated to a certain
//orientation (in degrees); length should be an odd number

void buildSE(pip::Image<unsigned char>& imSEUnRot, int orientation, int elementLength = 5){
        //First make the unrotated kernel
        pip::Dimension dim(elementLength,elementLength); // dimensions of SE
        imSEUnRot.resize(dim);
        imSEUnRot.fill(0);
        int lineHeight=elementLength/2;
        for(int i=0 ; i<elementLength ; ++i){
                imSEUnRot(i,lineHeight) = 1;
        }

        //Now rotate it over an angle orientation
        rotateNearestZ(imSEUnRot,orientation,(unsigned char)(0));
}

template <class Tin>
int RedLesionDetector<Tin>::run(const pip::Image<Tin>& imIn, pip::Image<Tin>& imOut){
        int errors = 0;
        // Do a closing with a linear structuring element at different
        //  rotational orientations and take the minimum of the resulting images.
        // Subtract the original from the result and put it in imOut (black top hat transform).
        vector<int> orientations;
        for(int i=0; i<nrOfOrientations; ++i)
        orientations.push_back((int)round(180.*i/nrOfOrientations));

        pip::Image<unsigned char> imSE;
        pip::Image<Tin> imClose;
        buildSE(imSE,orientations.at(0));
        grayClosing(imIn,imOut, imSE);
        for(int i=1 ; i<nrOfOrientations ; ++i){
                buildSE(imSE,orientations.at(i));
                grayClosing(imIn,imClose,imSE);
                for(int y=0; y<imIn.dim(1); ++y)
                for(int x=0; x<imIn.dim(0); ++x)
                imOut(x,y) = min(imOut(x,y),imClose(x,y));
        }
        for(int y=0; y<imIn.dim(1); ++y)
        for(int x=0; x<imIn.dim(0); ++x)
        imOut(x,y) -= imIn(x,y);

        return errors;
}

template<class Tin>
int RedLesionDetector<Tin>::run(const pip::Image<pip::Color<unsigned char> >& imIn, pip::Image<Tin>& imOut){
        pip::Image<Tin> imIn2;
        PreProcessFundus<Tin> ppf;
        ppf.setColorChannel('G');
        ppf.setMedianWidth(25);
        ppf.run(imIn,imIn2);

        RedLesionDetector<Tin> rld;
        rld.setNrOfOrientations(this->nrOfOrientations);
        rld.setElementLength(this->elementLength);
        int errors = rld.run(imIn2,imOut);
        return errors;
}
#endif